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Abstract 

Multistrain diseases have multiple distinct coexisting serotypes (strains). For some 
diseases, such as dengue fever, the serotypes interact by antibody-dependent en- 
hancement (ADE), in which infection with a single serotype is asymptomatic, but 
contact with a second serotype leads to higher viral load and greater infectivity. 
We present and analyze a dynamic compartmental model for multiple serotypes ex- 
hibiting ADE. Using center manifold techniques, we show how the dynamics rapidly 
collapses to a lower dimensional system. Using the constructed reduced model, we 
can explain previously observed synchrony between certain classes of primary and 
secondary infectives [14]. Additionally, we show numerically that the center manifold 
equations apply even to noisy systems. Both deterministic and stochastic versions of 
the model enable prediction of asymptomatic individuals that are difficult to track 
during an epidemic. We also show how this technique may be applicable to other 
multistrain disease models, such as those with cross-immunity. 
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1 Introduction 



In the study of infectious disease, a problem of interest is the coexistence 
of multiple strains. Examples of persistent co-circulating multistrain diseases 
include influenza [1], malaria [12], and dengue fever [11]. More recently, the 
avian flu viruses, including H5N1, were reported to coexist in several genotypes 
until 2004 [17]. Other examples of viruses possessing multiple strains may be 
seen in corona viruses, such as severe acute respiratory syndrome (SARS) [18]. 
The presence of multiple strains adds greater complication to models of disease 
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dynamics due to an increasing number of stages through different infection- 
recovery combinatorics. Reducing the dimensionality of the model is desirable, 
both to obtain a simpler model and to understand how the strains interact. 

In this paper, we focus on dengue fever, which has four co-circulating serotypes. 
It is believed that following infection with and recovery from one serotype, 
cross-reactive antibodies act to enhance the infectiousness of a subsequent 
infection by another serotype [16]. This phenomenon is termed antibody- 
dependent enhancement (ADE). Primary infections are less severe, often asymp- 
tomatic, but secondary infections are associated with more serious illness and 
greater risk for dengue hemorrhagic fever (DHF). 

Multistrain diseases with ADE-induced dynamics have been modeled previ- 
ously [10,14,5]. Models of multistrain steady state endemic behavior and its 
stability were considered for two strains in [9], where both mosquito popula- 
tion and humans were included. Although ADE was not explicitly modeled, 
the authors did consider different parameter values of contact susceptibility. 
Values of contact susceptibility less than unity model cross-immunity, whereas 
values greater than unity model increased susceptibility due to secondary in- 
fection. Ferguson et al. modeled two serotypes of dengue and showed that 
ADE can lead to oscillations and chaotic behavior [10]. Schwartz et al. devel- 
oped a similar model with non-overlapping compartments for greater ease of 
interpretation [14]. The bifurcation structure was obtained for both the au- 
tonomous model and a model with seasonal forcing. Chaotic solutions were 
found for a wide range of parameter values. It was observed from numerical 
simulations that outbreaks of the four serotypes could occur asynchronously 
but that certain primary and secondary infective compartments remained syn- 
chronized. In particular, all compartments currently infected with a particular 
strain exhibited in-phase outbreaks. Cummings et a/., using the same model, 
explored the competitive advantage that ADE confers to viruses. Other mod- 
els for dengue fever (e.g., [8]) have included the mosquito vector as in [9] but 
have not included the ADE effect. 

Several models of multistrain diseases with cross-immunity rather than ADE 
have appeared in the literature [1,6,8]. In these models, individuals who have 
recovered from infection with one serotype have reduced susceptibility to in- 
fection with other serotypes. Oscillatory dynamics have been observed. 

Summarizing the results to date, in models of secondary infections where ei- 
ther ADE or cross- immunity is considered, the steady state equilibrium is not 
always stable, and persistent oscillations may occur. For models with ADE, 
although ADE is the cause of oscillation onset, more complicated chaotic or 
stochastic behavior is observed over a much wider range of parameter val- 
ues than those of periodic oscillations [14]. Currently, there is no theory for 
predicting the interrelationship between primary and secondary infective com- 
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part merits. 

In this paper, we consider the problem of revealing the dynamical structure 
between the secondary and primary infectives. Specifically, we analyze the 
model presented [14] and demonstrate that the dynamics can be reduced to 
a lower dimensional model. Our analysis motivates the in-phase behavior pre- 
viously reported for primary and secondary infectives. The lower dimensional 
system allows the prediction of primary infective levels from knowledge of 
other compartment values. This result may be useful in disease monitoring 
because little epidemiological data exists for the asymptomatic primary infec- 
tives. We show numerically that our predictions hold even in noisy systems. 
We also discuss the relevance of this technique to other multistrain diseases 
that do not display ADE, such as influenza, in which infection with one strain 
yields partial immunity to other strains. 



2 Description of general n-serotype model 

We study the model of [14] for n co-circulating serotypes. An individual can 
contract each serotype only once and is assumed to gain immunity to all 
serotypes after infection from two distinct serotypes. This compartmental 
model divides the population into disjoint sets and follows the size of these 
compartments as a percentage of the whole population over time. The variable 
definitions are as follows: 

s Susceptible to all serotypes; 
Xi Primary infectious with serotype % 
Ti Primary recovered from serotype % 

Xij Secondary infectious, currently infected with serotype j, 
but previously had i (i ^ j). 

The model is a system of n 2 + n + 1 ordinary differential equations describing 
the rates of change of the population within each compartment: 





(2) 
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^1 = axi - pr-i J2\ x j + <A? Yl x kj ] - fJ'dU (3) 

' L> U - Pn I Xj + Oj x kj ) - a - /'./•''-./• (4) 



dt 



k+j 



The parameters fx, fi d , (3, and a denote birth, death, contact, and recovery 
rates, respectively. Antibody- dependent enhancement is governed by the pa- 
rameters (pi. Rates of infection due to primary infectives are of the form (3sXi, 
as in a standard SIR model, while infection rates due to secondary infectives 
are weighted by the appropriate 0, taking the form (3<f>jSXij. When = 1, there 
is no ADE, and both primary and secondary infectives are equally infectious. 
When > 1, the ADE appears as an enhancement factor in the nonlinear 
terms involving secondary infectives. 

Although the contact rate (3 could be time dependent (e.g., due to seasonal 
fluctuations in the mosquito vector population), we assume it to be constant 
here for simplicity. The fixed parameters throughout the paper are given by: 
fx = 0.02, fx d = or fx d = 0.02, (3 = 400, and a = 100, all with units of years -1 . 
These values are consistent with estimates used previously in modeling dengue 
[5,2]. We assume all serotypes have equal ADE factors: 0, = for all i. The 
ADE factor has not been measured from epidemiological data, so we test our 
results for various ADE values. 



The dynamics of Eq. 1-4 have been studied previously [14,2]. A bifurcation 
diagram in the ADE factor is given in Fig. 1 for the four serotype model. 
The dynamics are qualitatively similar for other values of n, the number of 
serotypes. For e [1,0 C ), the system has a stable endemic steady state. At 
a critical value, C , the system undergoes a Hopf bifurcation and begins 
to oscillate periodically. For slightly larger values, the periodic solutions 
become unstable and the system oscillates chaotically. Chaotic oscillations 
persist for most larger values for 0, with the exception of narrow windows of 
stable periodic solutions. 

It should be noted that previous studies of this model [14,2] used a death 
rate equal to the birth rate (na = fx = 0.02) so that the population remained 
constant in time. However, the endemic steady state cannot easily be obtained 
analytically when the mortality terms of Eqs. 1-4 are included. In our analysis, 
we omit the small mortality terms by setting ix d — 0. The dynamics with 
Hd = are qualitatively similar to that with Lid = 0.02 and have the same 
bifurcation structure. We demonstrate numerically that our results hold well 
even for the system with mortality. 

The larger parameters (3, a may be scaled by the small parameter fx, defining 
(3 = Po/fi, a = a /fx, (5) 
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Fig. 1. A bifurcation diagram for Eqs. 1-4 with n = 4, f3 = 400, a = 100, ji = 0.02, 
and Hd = 0.02. Shown is a Poincare section constructed from the extrema of a time 
series for susceptibles of length t = 100 years after a transient of t = 5000 years. The 
minima are denoted by grey points and the maxima are overlaid in black points. 
The Hopf bifurcation occurs at approximately (j) = 1.88333. 

where (3 , a are 0(1). Eqs. 1-4 with jj, d = have an endemic steady state of 
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3 Center manifold analysis 



Although the full n-serotype model has n 2 + n + 1 dimensions, the attractor 
lies in a much lower dimensional space. We claim that the attractor actually 
lies in In + 1 dimensions to a good approximation, due to a relationship 
between primary and secondary infectives. This result can be shown using 
center manifold analysis. We now outline a derivation for the case of two 
serotypes and state the center manifold equations for several other general 
cases of interest. 



3.1 n = 2 serotypes 



We begin with Eqs. 1-4 for two serotypes with no mortality (u^ = 0). We 
shift the variables so that the endemic fixed point (Eq. 6) is at the origin: 
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s = s — so, Xi = Xi — Xifi, Ti = Ti — r^o, Xij = Xij — Xijfl. The Jacobian in the 
shifted variables, evaluated at the origin, is given by 



1 + 
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to zeroth order in \i. 

The Jacobian is not diagonalizable, as it has only 5 linearly independent eigen- 
vectors. However, a new set of variables w is defined as follows: 



W = —J— 7 [X 21 - Xi, X 12 - X 2 , (1 + 0)S, X 1 + 0X 2 1, X 2 + 0Xi2, 
1 + 

0(X 21 - Xi) - (1 + 0)fi, 0(X12 - X 2 ) - (1 + 0)f 2 ] (8) 

The transformation matrix for this change of variables contains the 5 linearly 
independent eigenvectors of the Jacobian and 2 additional vectors selected to 
be linearly independent and put the system in the desired form. 

In the new variables, rescaling time to r = t//i, the dynamics are described 
by 



^p- = -a w 1 + P w 4 (w 7 + 0w 2 - w 3 ) + ^-^-(w 7 + 0u> 2 - w 3 ) (9) 

GST 2o"q 

dw 2 {3 

-a w 2 + f3 w 5 (w e + <pwi - w 3 ) + ^— -(w 6 + <jnvi - w 3 ) (10) 



(I: 2a 



-0 O (w 4 + w 5 ) - /3 (1 + 0)w 3 (w 4 + w 5 ) - ^-^-{\ + (j))w 3 (11) 
ctr 00 

G?U> 4 2 /i 2 /5 ( 



/5 O W 4 (0 2 W 2 + 0W 7 + W 3 ) + ^-^(0 2 W 2 + 0^7 + W3) (12) 

ctr 2cr 

-^0^5(0 ™1 + 0^6 + ^3) + "Ti 

20 O 

: cr (w 4 - iu 5 ) + /3 o u> 4 (-0 2 w 2 + 0w 3 - 0w 7 ) 



: /3 O W 5 (0 2 Wi + 0W 6 + W 3 ) + ^— ^(0 2 Wl + 0U> 6 + W 3 ) (13) 

ctr zao 

C?W6 2 



dr 
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+p w 5 [-0(1 + 4>)w 1 - (l + 4>)w Q ] 



2a 



W 2 + 0U>3 - (f)W 7 - 0(1 + 4>)wi - (1 + 0)w 6 



(14) 



dw 7 



■a {w 5 - w 4 ) + f3 w 4 [-0(1 + 4>)w 2 - (1 + (f>)w 7 ] 
+/3 w 5 (-(f) 2 w 1 + (pw 3 - (pw 6 ) 



2an 



-0(1 + 4>)W 2 - (1 + 0)U>7 - Wl + 0^3 - 0^6 



-i- = 



(15) 
(16) 



The Jacobian of Eqs. 9-15, evaluated at the origin and to zeroth order in //, is 
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(17) 



Eqs. 9-16 may be written in the compact form 



ii = Au + f (u, v, /i) 
v = Bv + g(u,v,/i) 
fi = 

where u = (wi,w 2 ), v = (w 3 , . . . ,wy), where A and B are constant matrices 
such that the eigenvalues of A have negative real parts and the eigenvalues 
of B have zero real parts, and f and g are second order in u, v, /i. Therefore, 
center manifold theory [3] can be applied. The system will rapidly collapse 
onto a lower dimensional manifold defined by 



w 1 = h(fi,w 3 ,w 4 , ...,w 7 ) 

W 2 = l(/l,W 3 ,W4, ...,W 7 ) 



(18) 
(19) 



We expand h and I to second order 
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7 7 7 7 

wi^h + ^2 hiWi + h^fi + Ki w i w 3 + h^Wi (20) 

i=3 i=3 j=i i=3 

7 7 7 7 

w 2 ~ /o + H ^ + ^ + EE + H W/^- ( 21 ) 

i=3 i=3 j=i i=3 

To solve for the coefficients of the expansion, we write equalities for and 
^p. This is done by equating the right hand sides of Eqs. 9-10 with the time 
derivatives of Eqs. 20-21, using Eqs. 11-15 to substitute for . . . , ^ and 
Eqs. 20-21 to substitute for wi, w 2 as needed. The coefficients of the expansion 
can be obtained by equating coefficients of like powers. 

The coefficients h , h iy h^, h^i and l , l iy l^, 1^ are all 0, while some of the h it j 
and lij are nonzero. After carrying out the above program, Eqs. 20-21 simplify 
to 



(22) 
(23) 

when the coefficients are substituted. 

Finally, we rewrite Eqs. 22-23 in the original variables. Thus we obtain the 
following approximation for the invariant manifold onto which the system 
collapses: 
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cr (xi - x 2 i) = A 

(7 (x2 - X 12 ) = A) 



2-0_ 30 

s ~- f2 + IT0* 2 + TT0* 12 

_ _ 2-0_ 30 _ 

s - ri + IT/ ;i + IT0 X21 



(Xi + 0X21 ) 

(x 2 + 0x 12 ) (24) 



The above equations for the center manifold are our main result of the paper. 
To bring out more of the structure, we make the following observations. We 
generally observed in numerical simulations that the infective compartments 
(and their deviations from the fixed point) are small compared to the suscep- 
tibles and recovereds. Thus the infective correction terms to s — fj may be 
dropped to obtain a simpler expression for the center manifold: 

cr (xi - x 2 i) = A (s ~ r 2 ) (xi + 0x 2 i) 

o-q (x 2 - x 12 ) = p (s - n)(x 2 + <j)x 12 ) (25) 

Given values for the susceptibles, primary recovereds, and secondary infec- 
tives, the primary infective compartments may be computed from the approx- 
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Fig. 2. Primary infectives with strain 1 vs. time for two-strain model with no mor- 
tality. Solid gray line: prediction from approximate center manifold equations (pre- 
dictions from Eqs. 24 and 25 lie on same curve); dashed black line: time trace from 
direct integration. The curves from Eqs. 24 and 25 coincide almost exactly. The dy- 
namics are chaotic with the parameter values used here: (3 = 400, a = 100, fi = 0.02, 
Hd = Q,<t> = 2-8. 

imate center manifold equations (either Eqs. 24 or 25). Figure 2 compares the 
center manifold prediction for primary infectives with that from simulations. 
(Results from Eqs. 24 and 25 are nearly indistinguishable.) The second order 
approximation to the center manifold leads to reasonable predictions for the 
primary infectives. In particular, we accurately predict the time and approxi- 
mate magnitude of the bursts, which correspond to epidemic outbreaks. The 
prediction for the infectives is occasionally negative because the center mani- 
fold equations are for the Xi, the deviations of infectives from their fixed point, 
and sometimes the predicted deviations are large enough that adding the fixed 
point Xifi yields negative values for the infective variables Xj. However, these 
errors are not physically important because they occur in quiescent regions 
where epidemic outbreaks are not occurring. Moreover, the timing of the pre- 
dicted outbreaks agrees well with actual outbreak time of the full system. 

3.1.1 Asymmetric ADE factors 

Antibody enhancement factors have not been measured experimentally, and it 
is possible that different serotypes could have different levels of enhancement. 
We address the possibility of unequal ADE factors for the two serotype case. 

A procedure similar to that used with symmetric ADE factors may be followed. 
In this case, an expansion for the location of the fixed point must also be used. 
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Fig. 3. Primary infectives with strain 1 vs. time for two-strain model with no mor- 
tality and asymmetric ADE factors. Solid gray line: prediction from approximate 
center manifold equations (Eq. 26 and 27); dashed black line: time trace from direct 
integration. Parameters: = 400, a = 100, u = 0.02, fi^ = 0, 0i = 2.5, 02 = 3.75. 
The ADE asymmetry measure is e = 0.5. 

For two strains with 0i = 0, 02 = 0(1 + e) arid |e| <C 1, the center manifold is 
approximated (to second order in the variables and /i, e) by 



x 1 - x 2 i + e (xi + 0x 2 i 
2 



r 2 + 



O"0 



1 + ' ' 1 + 

%2 - X12 - e (^2 + #12) 
2 



2(rr o + l)(l + 0)_ 

Xi + 0X 2 1 ) 



30 _ 

^2 + ~ 7X12 



(26) 



O"0 



2 (a + 1) (1 + 



= A) 



ri + 



-Xi + 



1+0 1 + 



£21 



(X 2 + 0X i2 ) 



(27) 



The center manifold equations again allow fairly accurate predictions of the 
primary infective compartments for small e (data not shown). Even when the 
|e| <S 1 assumption does not hold, the time and amplitude of bursts in the 
primary infectives are generally predicted accurately (Fig. 3). 



3.2 n > 2 strains 



It is important to extend the analysis to a larger number of serotypes. In par- 
ticular, there are four observed serotypes for dengue fever. The center manifold 
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technique may again be applied for n > 2 strains. In the expressions for the 
center manifold, it is convenient to define the sum of secondary infectives 
currently infected with strain k: 

n 

Zk = ( 28 ) 



For dengue fever, primary infectives are generally asymptomatic, and most 
hospital cases are secondary infectives. The current infecting serotype can 
be determined from serology measurements [13]. Zk, the sum of secondary 
infectives with serotype k, may be proportional to the number of hospital 
cases with serotype k. This quantity may be relevant to disease monitoring in 
addition to providing a shorthand notation for expressing the center manifold 
equations. 

We have completed center manifold analysis for n = 2,3,4 serotypes. The 
following equations for the center manifold summarize our results for n = 
2, 3, 4 and extrapolate them for general n: 



o"o W - z k ] = A) 
(To [(n - l)x jk - z k ]=p 
where 



i + k 



(n - l)rj - J2 r i + 9jk(x p , x pq , n, 0) 
i^k 



(x k + <pzk) (29) 
(x k + <t>z k ) (30) 



fk — 



9jk 



(n-l)(l + 0) \^ X: > ~ 



( 2n ~ W I V x-,-Tx, 



(31) 



n-(n-l)(f) 

(2n - 1)0 ; ^ _ ^ . 



(n-l)(l + 0) ^ ^ 



n 



(32) 



#3 



We observe that the terms fk and gjt may be neglected. They are linear func- 
tions of primary and secondary infective compartments. At the fixed point, the 
infective compartments are order /x, while the susceptibles and recovereds are 
order 1 (Eq. 6). It therefore is reasonable that deviations of the infectives from 
the fixed point, represented by the variables x p ,x pq , would generally be small 
compared to s, f\. This assumption generally holds true in our simulations. 
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For each of the n strains, Eqs. 29 and 30 provide n — 1 linearly independent 
equations, since some of the Eqs. 30 are linearly dependent, allowing n(n — 1) 
dimensions to be eliminated. Thus the dynamics of the system have dimension 
2n + 1. If the susceptible and primary recovered compartments are known, a 
single quantity z k (Eq. 28) for each strain is sufficient to describe the infective 
dynamics for that strain. Primary and secondary infectives can then easily 
be computed using Eqs. 29 and 30. We find that similar accuracy is attained 
when the f k , gj k are set to for convenience as with the f k , terms retained. 

Figure 4 compares the predictions from the center manifold equations (Eqs. 29 
and 30) with actual values for sample primary and secondary infectives. The 
full four-strain model with mortality (fid = 0.02) was used in the simulations. 
Predictions were made by assuming that susceptibles, primary recovereds, and 
sums of secondary infectives (z k ) were known. Although the derivation of the 
center manifold equations omitted mortality terms, they generate accurate 
predictions for the full model. 

3.2.1 Phase synchrony between compartments 

It has been observed in numerical simulations that primary and secondary in- 
fectives currently infected with the same serotype oscillate in phase synchrony 
[14]. Based on the center manifold reduction in dimension, we provide some 
analysis to motivate this effect. 

We convert the shifted variables (s, x iy f i^) to the original variables (s, x iy r iy x^) 
and use the center manifold equations 29-30 (with /, g terms omitted as ex- 
plained above) for substitution in the model ordinary differential equations, 
Eqs. 2 and 4. We obtain 



dt 



dx k 



az k + (3^2 r 'i ( x k + <pz k ) +)3 s-^rJ {x k ,o + <pz k ,o) (33) 




1) 



az k + (3^2n (x k + <pz k ) 




(34) 



where we have defined 



n 



z k = 



(35) 



i=l,ijtk 



for convenience. 



12 



(a) 




time (years) 



Fig. 4. (a) Primary infectives with strain 1 vs. time for four-strain model with 
mortality; (b) secondary infectives contracting strain 1 then 2 vs. time for four-strain 
model with mortality. Solid gray lines: prediction from approximate center manifold 
equations (Eqs. 29 and 30 with f k ,gj k terms neglected); dashed black lines: time 
traces from direct integration. The dynamics are chaotic with the parameter values 
used here: (3 = 400, a = 100, n = 0.02, /j, d = 0.02, r/> = 2.0. 

Taking the difference between Eqs. 33 and 34 and converting back to rescaled 
variables yields 

^f~( n ~ = I3[s-{n- l)fj] (x kfi + (j)z kfi ) (36) 

The difference between Eqs. 29 and 30 provides a useful substitution relation: 

a [x k - (n- l)x jk \ = (3 [s - (n - l)fj] (x k + (f)z k ) (37) 
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Substitution of Eq. 37 into 36 gives an equation for how the difference between 
primary and secondary infectives evolves in time: 



^0-(-.W01- ^y + y (38) 



where K = cr(x kt0 +(f)z k! o) is a positive constant of order /i. We restrict our con- 
sideration of Eq. 38 to those intervals on which x k (t) + (f>z k (t) is negative, thus 
avoiding singularities in [x k (t) + (j)Zk{t)\~ 1 . That is, we consider intervals in 
which the sum of primary and secondary infectives with serotype k (weighted 
by the ADE factor) is below the steady state value. In the chaotic regime, 
these intervals occur during dropouts when serotype k is present in low levels. 
Such intervals cover the majority of the time domain and are interrupted by 
outbreaks of serotype k. (Cf. Figure 4.) 



Integrating Eq. 38, we find that 



Xk(t) - (n - l)x jk {t) 

= [xk(to) - in - l)x jk (t )] exp < K J [x k (s) + ds > (39) 



From our previous assumption, the argument of the integral is negative and 
is of order 1/fJ 2 , since the infectives range between and the negative of 
the fixed point. Therefore, \x k (t) — (n — l)xj k (t)\ <C 1, and x k (t) and (n — 
l)xjk(t) approach each other rapidly, with their difference decreasing with 
exp [—0(1/ /j,)], where 1. 

We have demonstrated that primary and secondary infectives with serotype k 
are phase synchronized, differing only by an exponentially small term, during 
the intervals in which levels of serotype k are low. Outbreaks in the infectives 
occur in bursts, so although the relation \x k (t) — (n — l)xj k (t)\ <C 1 is not ex- 
pected to hold during outbreaks, the outbreaks have fast time scale dynamics. 
Therefore the duration of the outbreak is short enough that phase synchrony 
is not lost. It may be possible to examine the phase locking during the out- 
breaks, but that would require a detailed multiscale analysis, which is beyond 
the scope of the present paper. Even so, our results help to explain the phase 
synchrony previously observed between primary and secondary infectives cur- 
rently infected with the same strain [14]. 
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4 Stochastic perturbations 



We next consider Eqs. 1-4 with multiplicative noise and discuss whether 
Eqs. 29 and 30 can be applied to a noisy system. We include the noise in 
our numerical integration as follows. 

Let Y = {s,x i ,r i ,x i j} i j =l n and Y = F(Y). Define natural log coordinates 
y = {yi}, where ln(y) = y,i and y = f(y). Standard change of coordinates 
relates the two systems by: 

Y% — ^i/i( m Y), since | = y t (40) 



We add noise to the natural log coordinates: 

y = f(y)+77- (41) 



Here, 77 is a vector of random variables with mean zero and standard deviation 
a n . Transforming back to the original phase space, we obtain 



y i =y i / i (iBY) + 77 i y i (42) 
y = F,(Y) + ^y (43) 

Therefore, the noise is multiplicative and scaled by each component. Eq. 41 
was integrated numerically using a stochastic Euler method. 

The predictions from Eqs. 29 and 30 hold approximately true for a wide range 
of noise standard deviations a n . Sample results are shown in Figure 5. For the 
parameter values in Figure 5, the endemic steady state is stable for a n = 0. 
Adding noise perturbs the system away from the steady state, and as the 
noise standard deviation a n increases, the amplitude of the disease outbreaks 
increases. 

To quantify the accuracy of the center manifold prediction for various a n , we 
define a metric for the "distance" of a trajectory from the center manifold: 



d = ( I x i (* ) ~ x ^vred (t)\) , (44) 




where Xi tPre d is the predicted primary infective value using Eqs. 29 and 30 
(neglecting fk,gjk terms) and (.) t denotes a time average. We compute d by 
sampling every 0.01 year for 10 4 years (after first running for 10 4 years to 
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Fig. 5. Primary infectives with strain 1 vs. time for noisy four-strain model with no 
mortality, for various noise standard deviations. Solid gray lines: prediction from ap- 
proximate center manifold equations (Eqs. 29 and 30 with fk,gjk terms neglected); 
dashed black lines: time traces from direct integration, (a) a n = 0.01, (b) a n = 0.05, 
(c) a n = 0.1, (d) a n = 0.15. Other parameter values: (3 = 400, a = 100, /x = 0.02, 
^ = 0, 0=1.9. 



remove transients) in order to obtain good statistics. Sample results are given 
in Figure 6(a). The "distance" d from the center manifold, i.e., the error 
in the center manifold prediction, increases as the noise increases. However, 
the amplitude of the outbreaks also increases with larger noise, so it is not 
unexpected that d would increase. 

We additionally compute a normalized distance in lieu of a percent error so 
that the effect of noise may be more carefully assessed. (The actual percent 
error is not used because it would be greatly elevated by large percent errors 
in the dropout periods when the infectives are low.) The 10 4 year sampling 
period is divided into windows of 100 years, and the maximum value of the 
primary infective X\ (arbitrarily chosen) is computed for each window. These 
maxima give a scale for the outbreak amplitudes for each a n . The distance d 
is then divided by the average maximum for each a n to compute a normalized 
distance from the center manifold. Results are given in Figure 6(b). Using this 
metric, we find that the normalized error in the center manifold prediction 
actually decreases in noisy systems. 
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Fig. 6. (a) Average "distance" d from center manifold vs. noise standard deviation; 
(b) normalized distance vs. noise standard deviation, where distance and normalized 
d are as denned in the text. Parameters: (3 = 400, a = 100, \i = 0.02, = 0, <j> = 1.9. 

5 Conclusions and discussion 



We have analyzed a model for multistrain diseases with antibody-dependent 
enhancement. For a disease with n strains and two sequential infections, the 
model has n 2 +n+l equations. Using center manifold analysis, we have reduced 
the model to 2n + 1 dimensions. Although we have derived an approximation 
to the manifold on which the solutions for both deterministic and stochastic 
systems lie, if we use this approximation and evolve the dynamics of a model 
constrained to the manifold, we do not necessarily see the same bifurcation 
structure. This is due to the fact that the analysis itself is local and done near 
the steady state endemic point. Despite that limitation, the lower dimensional 
system has other uses which we have described here. 

We have partially explained the synchrony between primary and secondary 
infective compartments currently infected with the same serotype. During the 
intervals when a serotype is present in low levels, primary and secondary infec- 
tives of that serotype are approximately related by a constant of proportional- 
ity. Although our derivation does not hold during outbreaks of that serotype, 
the outbreaks occur in short bursts during which phase synchrony is not lost. 
As mentioned above, a full multiscale analytic treatment may lend more in- 
sight into the phase relationship between primary and secondary infections of 
a given serotype during the outbreak periods. 
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Our center manifold approximation enables prediction of each primary and 
secondary infective value if the disease parameters, susceptible and recovered 
compartments, and sums of secondary infectives (Eq. 28) are known. The pre- 
diction matches well with numerical simulations even when significant amounts 
of noise are added to the system. It has been observed [13] that most hospital 
cases of dengue fever are secondary infections. From serotype measurements of 
these hospital cases, estimation of the number of secondary infectives that have 
a particular serotype i should be possible. On the other hand, there is little 
epidemiological data on the frequency of primary infections because patients 
with primary infections tend not to be as seriously ill and are not admitted 
to hospitals. They are typically considered asymptomatic. Obtaining primary 
infective data for dengue monitoring purposes might require random sampling 
of large numbers of apparently healthy individuals, which would not be feasi- 
ble economically nor practically. Our center manifold equations might be used 
to improve monitoring of dengue outbreaks by predicting compartments that 
cannot readily be measured, such as the primary infectives. 

The fact that susceptible and recovered compartment data is necessary to 
make predictions for primary infectives is a potential drawback. However, we 
observe in numerical simulations that the susceptibles and recovereds vary 
more slowly than the infective compartments, which quickly spike during an 
outbreak. Therefore, it is possible that susceptible and recovered data could 
be obtained by sampling a population less frequently than would be required 
for infective data. If such sampling is possible, the center manifold equations 
will yield information about primary infectives from indirect measurements 
from population data. 

The applicability of our method to other models is also of interest. Multistrain 
models for diseases with cross-immunity, in which recovery from one strain 
leads to reduced susceptibility to further infection with other strains, have 
been developed [4,1,6,8]. We have analyzed a two-strain model with cross- 
immunity [15]: 




(45) 



(46) 



(47) 



(48) 
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where ip < 1 expresses the degree of cross-immunity and n — 2. This two-strain 
cross-immunity model is introduced in [4]. This model differs from the ADE 
model in that rather than giving the secondary infectives a different infectivity 
(weighted by 0), we give the primary recovereds a different susceptibility to 
further infection (weighted by ip) than the susceptible compartment. Using 
the coordinate transformation given in Eq. 8 (setting = 1), we can reduce 
Eqs. 45-48 from 7 to 5 equations and show a relationship between primary 
and secondary infectives. The eigenvalues of the Jacobian corresponding to 
rapidly collapsing directions are n{n — 1) eigenvalues equal to —a, the recov- 
ery rate. The center manifold analysis holds because the recovery rate is rapid 
compared to the birth rate. It seems reasonable to expect that other mod- 
els with cross-immunity might also be reduced and a relationship between 
primary, secondary, and perhaps higher order infectives might be derived. Be- 
cause cross-immunity models are often written in a different form with over- 
lapping compartments (e.g., [1,6,8]), some reformulation of the equations may 
be necessary before a center manifold approach can be used. 



The results presented here were primarily for a symmetric multistrain model 
with equal disease parameters for each strain. However, the case of unequal 
ADE factors was briefly considered, and a center manifold was found when 
the difference between ADE factors was a small parameter (Eqs. 26-27). To 
model a multistrain disease with realistic parameters, it may be necessary to 
relax the symmetry constraint and derive results for an asymmetric system. 
Such a derivation is unwieldy because the endemic steady state can no longer 
be determined analytically (except as an approximate expansion). New tech- 
niques may be necessary to extend our current approach to more asymmetric 
models. 



Another potentially interesting extension is to include spatial spread of the 
disease. It has been observed from epidemiological data in Thailand that while 
dengue outbreaks within a single town may be of primarily one serotype, con- 
current outbreaks in another town may have a different serotype predominant 
[7]. A similar center manifold analysis might be attempted on a patch model, 
in which the population is divided into several spatially distinct patches that 
are coupled. 
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